function [cd, rho]=CD(e)
% e is NxT
[N,Tmax]=size(e);
pom=0;


cd=0;   % CD test without averages
rhos=0;
for i=1:N-1
    for j=i+1:N
        e1=e(i,:)';
        e2=e(j,:)';
        ind=not(isnan(e1.*e2));
        Tij=sum(ind);
        ui=e1(ind); uib=sum(ui)/Tij; uid=ui-uib;
        uj=e2(ind); ujb=sum(uj)/Tij; ujd=uj-ujb;
        vi=(uid'*uid)^0.5;
        vj=(ujd'*ujd)^0.5;
        dij=uid'*ujd;
        rij=dij/(vi*vj);
        cd=cd+((Tij)^0.5)*rij;
        rhos=rhos+rij;
        pom=pom+1;
    end %j
end %i
rho=rhos/pom;
cd=cd*( (2/N/(N-1))^0.5);
end

